home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / Libraries / Advanced I⁄O v2.3 / Advanced i⁄o / histogram.cc < prev    next >
Text File  |  1995-05-29  |  4KB  |  135 lines
  1. // This may look like C code, but it is really -*- C++ -*-
  2. /*
  3.  ************************************************************************
  4.  *
  5.  *            Computing the Histogram  
  6.  *            and related operations
  7.  * $Id$
  8.  *
  9.  ************************************************************************
  10.  */
  11.  
  12. #ifdef __GNUC__
  13. #pragma implementation
  14. #endif
  15.  
  16. #include "histogram.h"
  17. #include "myenv.h"
  18. #include "std.h"
  19. #include <math.h>
  20. #ifndef M_LOG2E
  21. #define M_LOG2E         1.4426950408889634074
  22. #endif
  23.  
  24.                 // Constructor: constructs a clean histogram
  25.                 // ready to accumulate statistics for symbols
  26.                 // within [lwb,upb]
  27. Histogram::Histogram(const Histogram::Symbol lwb, const Histogram::Symbol upb)
  28.      : symbol_lwb(lwb), symbol_upb(upb),
  29.        no_potential_symbols(upb-lwb+1)
  30. {
  31.   assure(no_potential_symbols > 1,
  32.      "Bracketing interval for input symbols should include at least "
  33.      "two symbols");
  34.   counts = new unsigned int[no_potential_symbols];
  35.   memset(counts,0,sizeof(counts[0])*no_potential_symbols);
  36. }
  37.  
  38.                 // Destructor
  39. Histogram::~Histogram(void)
  40. {
  41.   assert( counts != 0 );
  42.   delete [] counts;
  43. }
  44.  
  45.                 // Count a symbol
  46. void Histogram::put(const Symbol symbol)
  47. {
  48.   if( symbol < symbol_lwb || symbol > symbol_upb )
  49.     _error("Symbol %d is out of interval [%d,%d]",
  50.        symbol, symbol_lwb, symbol_upb);
  51.   counts[symbol-symbol_lwb] += 1;
  52. }
  53.  
  54.                 // Count a symbol coercing it to the expected
  55.                 // range if necessary
  56. void Histogram::put_coerce(const Symbol symbol)
  57. {
  58.   counts[ ( symbol < symbol_lwb ? symbol_lwb :
  59.           symbol > symbol_upb ? symbol_upb :
  60.           symbol ) - symbol_lwb ] += 1;
  61. }
  62.  
  63.                 // Return the count for a symbol
  64. int Histogram::get(const Symbol symbol) const
  65. {
  66.   if( symbol < symbol_lwb || symbol > symbol_upb )
  67.     _error("Symbol %d is out of interval [%d,%d]",
  68.        symbol, symbol_lwb, symbol_upb);
  69.   return counts[symbol-symbol_lwb];
  70. }
  71.  
  72.                 // Give the no. of distinct symbols, ie
  73.                 // symbols with non-zero count
  74. int Histogram::no_distinct_symbols(void) const
  75. {
  76.   register int no = 0;
  77.   register unsigned int *cp;
  78.   for(cp=&counts[0]; cp<&counts[no_potential_symbols]; )
  79.     if( *cp++ != 0 )
  80.       no++;
  81.   return no;
  82. }
  83.  
  84.                     // Print the histogram
  85. void Histogram::print(const char * title) const
  86. {
  87.   message("\nHistogram '%s'\n"
  88.       "\nValue    Quantity\n",title);
  89.   register int i;
  90.   register int no_symbols = 0;
  91.   register int total = 0;
  92.   register unsigned int *cp;
  93.   for(i=symbol_lwb, cp=&counts[0]; i <= symbol_upb; i++,cp++)
  94.     if( *cp != 0 )
  95.     {
  96.       total += *cp;
  97.       no_symbols++;
  98.       message("%2d     %5d\n",i,*cp);
  99.     }
  100.   assert( cp == &counts[no_potential_symbols] );
  101.   message("\n****** Total no. of symbols %d,  total count %d\n\n",
  102.       no_symbols,total);
  103. }
  104.  
  105.  
  106.             // Estimate a zero-order entropy of a source
  107.             // represented by the histogram
  108.             // Using Shannon formula
  109.             // log2(f) - SUM[ fi*log2(fi) ]/f
  110.             // where fi is a count for symbol i, f is a total
  111.             // count
  112.             // The formula tells how many bits are necessary
  113.             // to represent a single *average* symbol emitted
  114.             // by the (stationary) source
  115.             
  116. static inline double log2(const double x)  { return log(x) * M_LOG2E; }
  117.  
  118. double Histogram::estimate_entropy(void) const
  119. {
  120.   register long total_count = 0;
  121.   register unsigned int *cp;
  122.   for(cp=&counts[0]; cp<&counts[no_potential_symbols]; cp++)
  123.     if( *cp != 0 )
  124.       total_count += *cp;
  125.   
  126.   assert( total_count > 0 );
  127.   
  128.   double accum_log = 0;
  129.   for(cp=&counts[0]; cp<&counts[no_potential_symbols]; cp++)
  130.     if( *cp != 0 )
  131.       accum_log += *cp * log2( *cp );
  132.  
  133.   return log2(total_count) - accum_log/total_count;          
  134. }
  135.